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00 ■ Abstract 

<n : 

Considering the Gross-Pitaevskii integral equation we are able to formally obtain an analytical 
solution for the order parameter $(#) and for the chemical potential \i as a function of a unique 

: 

dimensionless non-linear parameter A. We report solutions for different range of values for the 
repulsive and the attractive non-linear interactions in the condensate. Also, we study a bright 
soliton-like variational solution for the order parameter for positive and negative values of A. 
Introducing an accumulated error function we have performed a quantitative analysis with other 

o : 

O . well-established methods as: the perturbation theory, the Thomas-Fermi approximation, and the 

\ numerical solution. This study gives a very useful result establishing the universal range of the 

' A-values where each solution can be easily implemented. In particular we showed that for A < —9, 

-xt- : 

the bright soliton function reproduces the exact solution of GPE wave function. 
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I. INTRODUCTION 



Since the unambiguous experimental realization in dilute ultra-cold atom cloud of the 
Bose-Einstein condensed phase (BEC)^>2jM^&Z a lot of work has been devoted for search- 
ing the dynamic and physical properties of the nonlinear matter waves and excitations of 
the condensate. The observation of this effect in dilute atomic gases has allowed to invoke 
the weakly interacting mean field theory to describe the properties of the BEC systems.- 1 ^ 
Hence, the dynamics of the process for the order parameter has been ruled by equations 
of nonlinear Schrodinger (NLS) type and mainly by the Gross-Pitaevskii equation (GPE). 
Nowadays, NLS equations with attractive (negative scattering length}^ and repulsive (pos- 
itive scattering length)Ai nonlinear interactions have been reported to describe experimental 
observations of different types of wave solitons.— Most of the theoretical work has been 
devoted to implement numerical solutions of the GPE for the order parameter (see Ref. 9| 
and references therein). To study and to control the physical properties of the condensate 
it will be very useful to manipulate analytical expressions for the chemical potential and 
for the order parameter as well. A typical example of the great physical interest is the 
attention devoted to the collective excitation spectrum of a BEC. In this case we have to 
deal with the time-dependent GP equation under the linear response approximation. Here, 
as input parameter, we have to insert in the Bogoliubov equations^ the order parameter 
and the chemical potential solution of the NLS. In that sense, some analytical results for 
those magnitudes are priceless. Also, many order problems can be stressed on the field of 
cold atoms BEC as the dynamical stability— atomic current in an optical lattice,— etc. 

In order to achieve closed solution the variational procedure with a Gaussian as a trial 
wave function has been proposed (see Ref. 16| and reference therein). Nevertheless, it is 
well know that this Ansatz does not reproduce well the properties of the condensate. For 
example, in the repulsive interaction case and in the strong nonlinear limit, the shape of the 
order parameter should be similar to the Thomas- Fermi (TF) solution (see below Eq. (j4j)). 
Moreover, the results obtained with the standard variational procedure is, in many cases, 
qualitative and, even if the real shape of the wave function resembles the trial wave function, 
the variational method is not always a good reference for solving nonlinear equations.— 

Nowadays, available numerical methods for solving differential equations are fast and 
accurate. Nevertheless, if the evaluation of several physical magnitudes is carried out, such 
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as the optical properties among others, or to control the properties of the condensate (as we 
just discussed above), this advantage is lost due to cumbersome numerical computational 
procedures that must be performed at the end of the calculation. Moreover, if we work with 
a given basis of functions, it is difficult to know a priori if, in fact, the basis is a complete 
set for the Hilbert space of the specific nonlinear equation. Also, the type as well as the 
swiftness of convergence to the real solution is not always well established. In that sense, to 
implement manageable analytical expressions for the order parameter, where the accuracy 
and the absolute error of the obtained solution are controlled, has becomes a necessity. This 
is a fact in the study of nonlinear equation and in particular for the GPE. 

A description of the order parameter in terms of a controlled truncated basis becomes 
a useful tool if we are dealing with not many implemented functions and the degree of 
accuracy is well established. So, the obtained expansion will be given by a sum of few basic 
functions, allowing in that way to handle with explicit solution to describe the physical 
properties of the condensate. Unfortunately, the beauty of such a mathematical result is 
restricted to certain range of values for the parameter involved in the nonlinear equation 
under study. The challenge is to find precisely this range of convergence, to give the absolute 
error in terms of physical parameters, and to provide other handled compact solutions outside 
the obtained range of the desire accuracy. We would like to remark that the most important 
requirements for analytical solutions are simplicity, flexibility, and the viability to be used 
in perturbation approaches for the calculations of physical properties. 

In this paper we present different methods of solutions of the time independent GPE based 
on the equivalent integral GPE and its relation with the Green function of the corresponding 
linear operator, on the soliton solution, and on a bright soliton-like variational function. This 
discussion will provide general analytical expressions for the order parameter and for the 
chemical potential in a universal range of the non-linear interaction parameter. 

To describe the order parameter we started with the isomorphic one-dimensional 
nonlinear Gross-Pitaevskii equation (GPE), 8 ' 9 which can be written as 

— + -mw 2 r$ + A $ $ = ^$ 1 

2m dx l 2 

with the normalization condition 
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dx |$| 



(2) 



In the above equation \x represents the chemical potential, u is the trap oscillator frequency, 
m is the alkaline atom mass, and A is a self-interaction parameter describing the interaction 
between the particles. 

Equation ([1]) presents an explicit solution if the non-linear term (\ |$| 2 ) is larger than 
the mean value of kinetic energy operator. This approximation, known as TF,— ^ provides 
simple expressions for the chemical potential and the wave function given by 
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(4) 



where A = X/l hu>, l = ^Jh/moj, and the value of A > is restricted by Eq. (jlj. 

The next section is devoted to develop the proposed methods to solve Eq. ([I]), beyond 
the above typical TF approximation, . The main goal is to obtain explicit representations 
for the whole range of the self-interaction parameter A (negative and positive values) and to 
show the range of validity for each particular method of solution. 



II. ANALYTICAL APPROACHES 



First we will study the variational method based on a soliton wave function as Ansatz 
function, secondly we analyze the validity of the spectral method based on the equivalency 
between the integral and differential equation ([TJ and the Green function, solution of the 
linear harmonic oscillator operator. Moreover, using the obtained general formalism we 
report perturbation solutions for $ and fi in terms of the non-linear parameter A. For 
sake of comparison and in order to check the accuracy of the implemented approaches, the 
numerical solution of Eq. ([!]) is also addressed. 



A. Variational method: Soliton approach 



The variational method, valid for positive as well as negative values of A, could provide a 
simple picture of the main physical characteristics of the BEC. Without the trap potential, 
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the GPE ([I]) reduces to the nonlinear Schrodinger equation which for A < admits the 
stationary normalized bright soliton solution 

$s(x,K)= (-) sech(Kx). (5) 
Here, the chemical potential fis and the inverse of the soliton length K are expressed by 

m\ 2 m\ 

In order to solve Eq. ([1]) for all values of A we propose as variational Ansatz the bright 
soliton (jSJ) where K is taken as a variational parameter. The Ritz's variational method 
applied to the NLS ([I]) provides for the chemical potential fi(K) the parametric equation 
(see Appendix A) 

h 2 K 2 mu 2 K\ 
Vvar(K) = -^-a + ^(3 + — 7 , (7) 

where K must fulfill the dimensionless equation 



b 4 + b 3 - 5 = (8) 

with b = Kh 2 /(m\) and 5 = (l/A) 4 7r 2 /4. Accordingly, Eq. (JTj) is reduced to the simple 
relation 

tZL = _ ^) (9) 

nijj 6 b z 

and for the order parameter we get 

r- /A6\ 1/2 x 

Vk^var = I -J J Sech(Ab-). (10) 

Equation (jSJ) is a fourth-degree algebraic equation with only one real physical meaningful 
solution, which depends on the sign of the non-linear interaction parameter A. In order to 
get a more clear view of the solution for Eq. (jSJ), we carry out separate calculations at A = 
and for the strong repulsive (attractive) limit A — >• oo (A — > — oo). 
If A -> we obtain from Eqs. (jEJ) and (JUj) 

6xA=J| (11) 
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FIG. 1: Variational parameter b as a function of A. The asymptotic limits (I13p and (|15p . and the 
behavior at A ~ 0, Eq. (jlip . are indicated by dashed lines. 



and 

/W(A=0) 7T 



(12) 



fojj 6 

In the strong attractive limit (A <C — 1) and keeping the leading term in Eq. ([8]), the possible 
physical solution has the asymptotic behavior 

6 (A - -oo) = 1 + o - y (13) 

and the chemical potential \i is given by 

/w(A -> -oo) 1 2 



ftc<j 6 
In the repulsive limit case, A ^> 1, Eq. (jSJ) yields 



— . (14) 
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with 

(A — > OO) / 7T\/2 



(A) 1 - (16) 
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We can compare the above limit solutions with those obtained by the TF approximation, 
Eq. (j3J), and the exact soliton solution, Eq. (jSJ). The relative errors are equal to 



oo 



Htf(oo) 



2 

7T\ 3 



1 » 0.0312. 



(17) 



and 



-oo 



/xs(oo) 



yUS^OOJ 



0.3333. 



(18) 



From the above relations we conclude that the variational wave function (jSJ) provides a 
better solution of the GPE in the strong repulsive case than for the attractive one. At A = 
a relative error of 0.0472 is reached by comparing Eq. (112j) with the exact solution of the 
harmonic oscillator problem /i/hu = 0.5. In Fig. 1 we present the parameter b, solution of 
the Eq. flS}, as a function of the dimensionless parameter A. Also, the limit solutions are 
indicated by dashed lines. It can be seen that the calculated asymptotic behaviors (fTTj) . ({TBI) , 
and (|T5|) at A = 0, A — > — oo, and A — > oo, respectively, are quickly reached by the exact 
solutions of Eq. (jSJ). It is not surprising that the Ritz's variational method failed to get a 
closed analytical solution of the differential GPE. The variational method here implemented 
is only valid for linear differential equations or the corresponding Lagrangian of the problem. 



B. GP integral equation: Green function solution 

One of the most powerful analytical method used to solve differential and integral equa- 
tions corresponds to the Green function formalism (GFF). In order to implement this math- 
ematical technique to the non-linear Schrodinger equation we rewrite ([T|) as 

M*l = -££ + 5™-V* = /<*)■ P») 

Here f(x) will be considered as an inhomogeneity in the differential equation and equal to 

f(x) = ( f i-xmx)\ 2 Mx). 

Function $(x), solution of Eq. (|T9|) . can be cast in terms of the Green function G(x,x f ) of 
the linear operator L [<&] . Formally, we can write $( function of the inhomogeneity 

f{x) as 20 
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/oo 
G(x,x')(fi- \\$(x')\ 2 )<5>(x')dx'. (20) 
-oo 

The above expression corresponds to the GP integral equation for the order parameter $(x) 
We observe that the integral equation ( 120]) has a symmetric kernel, G(x,x'), which fulfills 
the differential equation 

Lq [G(x, x')} — 5(x — x') . 

To write the formal solution ( 1201) in terms of the Green function of the operator Lq, the 
function f(x) has some constrains.— 1 ^ In our case, all functions and the Green function 
also, have to fulfill the boundary condition $(x) — > as x — > ±oo.. This guarantees that 
the inhomogeneity f(x) belongs to the same Hilbert space of the linear operator Lq. The 
kernel G(x, x') is the given by the following spectral representation 



ra=0 



Huj(n+ 1/2) 

with <p n {x) being the harmonic oscillator wave function^ 



] \ ( —x 2 \ „ / x 



^ = {^¥m ) ^{w) H "\tJ- (22) 

We have to note that according to the general theory of Fredholm integral equations,— 1 ^ 
the set of functions appearing in the spectral representation of a symmetric kernel, {ip n (x)} 
in the present case, represents a complete set of functions for the given Hilbert space of the 
GP integral equation (T2"U|) . Hence, the convergence of the expansion f[2Tl) is guaranteed and 
we can insert the spectral representation of G(x, x') in (TSUI) and interchange the integral and 
infinity expansion (l2Tj) . Thus 

$ = E tJ[n + \/2) I Mx ' )ifI ~ X \^ x ')\ 2 )^ x ') dx ' ■ ( 23 ) 
From ( 1231) it is straightforward that the general solution for the order parameter $ has an 
explicit representation through the harmonic oscillator if n {x) as 

oo 

j> n (x)C n (Ai). (24) 

n=0 

Since the inhomogeneity f(x) belongs to the same Hilbert space of the symmetric kernel of 
the of Fredholm integral equation (120]) . the convergency of the series ( )24|) in energy to the 
function $ is guaranteed.— 
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In the present case the coefficients C n (/x) are restricted to obey the relation 

C n = J fa^+yzf nWO* - A mx'tfMx'W. (25) 

Inserting the convergent series (1241) in Eq. ( 1251) . it follows that the vector coefficient C(/x) 
must fulfil the non-linear equation system 

[A(ji) + ACT ■ C] C = 0, (26) 

where 

™ + ~-^W- (27) 



2 

and T p i mn is a fourth dimensional matrix defined in the Appendix B. 

The order parameter $ and the chemical potential as a function of the dimensionless 
parameter A are obtained by solving the non-linear equation system f l26|) . Although the 
mathematical complexity of Eq. ([TJ) has been reduced, Eq. (|26|) is nevertheless an infinite 
generalized eigenvalue problem for //(A) and C(/i(A)). The complexity of the problem 
depends on the sign and the values of the non-lineal parameter A but the key issue is how 
quickly converges the series in (|24|) or equivalently, the non-linear equation system (|26[) . 
This important problem is addressed in the next section. 

We have to mention that the obtained problem (T26]) is isomorphic to Galerkin method. 
The former one is a generalized variational method where for a given equation L[F] = 
L [F] + L p [F] it is possible to choose a certain basis {g^} of the operator L and to expand 
the function F in term of the given basis. The choice of the operator L (which must 
include the boundary conditions) is not unique and certain degree of freedom prevails. To 
guarantee that the expansion converge to the real solution, the picked out operators L and 



L p have to fulfil certain mathematical conditions (see Ref. 



24j for a detailed description of 



this mathematical treatment). This crucial question is not trivial when we are dealing with 
non-linear equations as the NLS. In our case the mathematical treatment above developed is 
based on the properties of the Fredholm integral equations and can be considered a rigorous 
demonstration of the validity of the expansion (1241) and the convergence to the correct 
solution. 
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We obtained the ground state solution <3> in terms of a truncated basis set, {(p n (x)} (n = 
1, .../), by defining the finite dimensional nonlinear Hill determinant eigenvalue equation^ 

||A<( / )( / i(A),C)||=0, (28) 

where A^nm, (n,m = 1,2, ...I) are the corresponding matrix elements according to the Eq. 
(126]) . Since the scaling of any direct numerical algorithm of integration implemented to 
obtain the tensor T p i mn is of the other J 4 x P (P is the number of grid points) the numerical 
implementation becomes a cumbersome task and non-efficient method of evaluation. To 
get a better efficient algorithm than those based on a direct numerical integration of the 
tensor T p i mn , it is necessary to exploit its analytical representation together with its sym- 
metry properties. This analysis is presented in the Appendix B allowing a straightforward 
evaluation of the tensor T p i mn . 

To solve Eq. (128]) . we have implemented the Neumann iterative procedure in a finite 
basis of dimension /. For a given iteration and since the functions {ip n (x)} define a complete 
set for the GPE, obeying the natural boundary conditions, <p n (x) — > for x— > ±oo, the 
roots of the determinant (T281 converge to the exact ground state solution of Eq. (EQ)) 
and lim/^oo //^(A) = /i(A). 25 The numerical procedure starts from a trial vector C ^ C 
and iteratively we obtain the k — th approximation. In each step, the matrix (128]) must 
be recalculated by using the new eigenvector C ^ C. The procedure is repeated until 
< 5 C and (or alternatively) l/i® — < fko ■ 5^, where 5 C and 5 M are 

the desirable accuracies for the coefficients and the chemical potential, respectively. For the 
iterative procedure, it is useful to introduce a control parameter e G [0, 1], so that 

C n =yJe(Ck k - 1) ) 2 + (l-e)(Ck k) ) 2 . 

This procedure is faster and accurate for positive and small negative values of the non-linear 
parameter A > —5. For A < —5 however, the size of the matrix we have to deal with grows 
as |A| does. In the former case, a basis set of 25 functions allows at least 5 significant figures 
in the calculation of //, while for the later at least 50 oscillator wave functions {<f n } were 
sorted in order to reach the same accuracy at A = —10. 

According to the Neumann iterative procedure we have to introduce an initial starting 
q(o) vector. This vector can be chosen according to the desirable A value and the following 



r (k) _ r (k-i) 
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criteria can be established: i) For the dimensionless interaction parameters |A| < 1.5, the 
coefficients Cn)m = 5 n ^ m . ii) If A > 5 the asymptotic limit of the TF approximation wave 
function given by (jlj) is a good starting iterative procedure, iii) For attractive interaction 
and A < — 1.5 the soliton wave function approach (jSJ) is useful as initial condition . 



C. Perturbation theory 

It is useful to get expressions for [i and the order parameter $ through a perturbation 
approach since these are easily handled solutions. Also, the explicit perturbation expressions 
can be implemented as a method to control other solutions in particular the numerical ones. 
If the nonlinear term H p = A |$| is considered as a perturbation in comparison to the trap 
potential mwV/2, the chemical potential and the vector C in Eq. (1261) can be sought in 
the form of series, i.e. 

c m = c<» + \c£> + \*c<» + ...., 

ii = ^ + A M « + A 2 /i (2) + 



Taking only the second order interaction in A, Eq. ( l26l) yields 



„ = £ + £ Tra _ 3 gy!!W (29) 



Using the properties of the matrix T p i mn given in the Appendix B we get 

fi 1 A 3 a2 <A (2m- 1) 



J^ A 2 V- ^<'t-,y. 

9^ A 2^ 04m(rr,\\2 ■ ^ 



huj 2 2tt ^ 2 4m (m!) 2 

v m=l v ' 



Using that 



(2m)! -91 (i 2x2 + 1 I 

2 3m (m!) 2 m(x 2 + l) m ~ ~ 11 I 2y 2(x 2 + 1) + 2 

the chemical potential up to second order is reduced to the following useful expression 

= - + -jL - 0.033106 x A 2 . (31) 



Finally, the normalized order parameter $ including terms to the first order can be expressed 
as 

v m=l v ' 
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FIG. 2: (Color online) Normalized order parameter \/Tq&(x/Iq) for the positive dimensionless self- 
interaction A values: a) 25, b) 15, c) 10, and d) 5. Solid line: Solution (|24p . Dashed line: Soliton 
variational approach. Dot: Thomas-Fermi approximation. Empty circles: Numerical solution. 

D. Numerical Solution 

A comparison of the obtained analytical solutions with direct numerical calculations is 
an important control for validating the mathematical methods here introduced. In order to 
solve ([1]) numerically, we choose a finite difference method where for the second derivative 
we select a simple three-points approximation with uniform spacing, so that the differential 
equation can be rewritten as a symmetrical tri-diagonal matrix. The eigenvalue problem for 
the obtained matrix can then be solved by the usual methods. Explicitly we have 
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FIG. 3: (Color online) Accumulated error function rj for the repulsive interaction as a function of 
A. Solution (|24p (solid line), soliton variational approach (|10p (dashed line), and Thomas- Fermi 
function @ (dot line). Inset: Perturbation wave function (|32p (dot-dashed line). 
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where 
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2/ 

-r 2 , 



(34) 



$( — 2^ + (i — 1)5), i = 1,2, ...L, and 5 is the discreet step. The presence 



of the wave function inside the matrix in the left side of Eq. fl33l) enforces the use of some 
kind of iteration procedure in order to solve the non-linear problem. That is, for the vector 
v of components f« (see Eq. (13~4"|) ) we set 



0,1,2,..., 



v [F<*>] , k 

where is a certain trial function. We started with certain F^(x) 



(35) 



$(°)(x) evaluated 

at the Xi mesh points. After that, we find the approximate eigenvector $ and eigenvalue 
/i of the ground state solution of Eq. (1331) . The new trial function is obtained by the 
expression 



(fc-i) 



'1-e) 



with e G [0, 1]. This procedure is repeated until 



< 5$ and (or alternatively) 



[/i^ — /i^ -1 -*! < /icj • 5^, where 5$ and <5 /t are the desirable accuracies for the wave function 
and t ie chemical potential, respectively. A similar procedure has been used with success in 
Ref. 26( for a two component BEC. The practical implementation of the above described 
method is mainly straightforward, however, due to the influence of the non-linear term, 
the accuracy and speed of convergence is critically dependent on the correct choice of the 
parameter e. Our experience shows that the best value e depends on the value and sign of 
the non-linear term. 



III. RESULTS 

We shall now discuss the accuracy and the reliability of the above implemented methods 
of solution, by studying independently the repulsive and the attractive interaction cases. 
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FIG. 4: (Color online) Chemical potential in units of huj as a function of dimensionless self- 
interaction parameter A. Solid line: Equation (|26l) . Dashed line: Soliton variational approach 
([9]). Dotted line: Thomas-Fermi approximation ([3]). Empty dot: Numerical solution. Inset: 
Perturbation theory (|3ip (dot-dashed line). 

A. Repulsive interaction 

Figure 2 displays the order parameter \flo§(x / Iq) for several values of A. The variational 
solution of the Eqs. ([9]) and fTTUT) is represented by dashed lines, solid lines present the 
calculation using Eq. (I2"4l . and the TF approach, following Eq. (jl]), is indicated by dots. 
Empty circles show the obtained numerical solutions of the GPE ([T|) following the procedure 
described in Sec. II D. It can be seen that the wave function is more delocalized and the 
maximum of Q(x/Iq) decreases as A increases, thus the condensate spreads as the non- 
linear term increases. Also, in Fig. 2, the differences between all calculated analytical 
representations and the numerical procedure are qualitatively displayed. As already known, 
the TF approach reproduces well the properties of the condensate for large values of A, while 
the proposed variational solution exhibits a better approximation for small values of A. In 
general, we have obtained very good agreement between the solution fl24|) and the numerical 
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FIG. 5: (Color online) The same as Fig. [2] for the attractive dimensionless self-interaction A values: 
a) -20, b) -15, c) -10, and d) -5. Dotted line represents the soliton solution (|5|). 

solution for all considered values of dimensionless interaction parameter A. Nevertheless, 
it is useful to define a magnitude that quantify the quality of the implemented analytical 
solutions. Hence, we have introduced the accumulated error function 

/oo 
\®num(x) ~ $i(aO| dx, (36) 
-co 

where <& num is the numerical solution of Eq. (1). The above magnitude gives a direct 
estimation of the total error introduced throughout the whole interval — oo < x < oo. Since 
in each given point x G (— oo, oo) we add the modulus of the difference between <& num {x) and 
$i(x), then rji determines the maximum accumulated error for the analytical wave function 
$i(x). Figure 3 presents the estimated error rji as a function of the dimensionless interaction 
term A for all functions considered: the Thomas- Fermi (dot line), the solution ff24l) (solid 
line), and the soliton variational approach (dash line). From the figure it can be seen that 
the best analytical solution (77 < 0.033) is reached by using ff2M . while the TF approximation 
approaches, asymptotically, to the exact solution . The soliton variational solution exhibits 
its minimum error (rj < 0.2) for A < 2, reaching a maximum error at A « 10. One can notice 
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FIG. 6: (Color online) The sane as Fig. [3] for the attractive interaction as a function of A. Dotted 
line represents the soliton solution ([5]). Inset: Accumulated error function rj for the perturbation 
wave function (|32p (dot-dashed line). 

that $„ ar is a better approach than the TF for A < 3.6. The accumulative error introduced 
by the perturbation wave function (dash-dot line) is also shown in the inset. In general, the 
accuracy of the series (|24p can be greatly improved if large matrixes are implemented. In 
our calculations a few functions (a 50 x 50 matrix) was necessary to achieve an accuracy of 
10~ 8 for the chemical potential. In the case of the numerical procedure (133]) . values of fi/fko 
were calculated with an uncertainty of 10 -10 . 

In Fig. 4, we compare the calculated chemical potential in units of the energy trap Ttw 
according to the analytical methods outlined in the previous section. The Thomas-Fermi 
approach following Eq. (jSJ) is indicated by dots, the soliton variational solution obtained by 
solving the Eqs. dHJ and (Q is represented by a dashed line, while the solid line presents 
the calculation using the Hill determinant fl28|) . In the inset, the comparison with the 
perturbation theory given by Eq. fl31~|) (dash dot line) is also shown. The numerical solution 
is also presented by empty dots. As expected the TF limit increases its accuracy, i.e., less 
than 3% of error at A = 10, as the non-linear parameter increases. No differences can be 
observed in the scale of the figure between the numerical solution and the chemical potential 
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using Eq. f|28|) . The soliton variational calculation presents a larger error for A > 7. Also, 
in the figure we can observe the relative error of 0.0312 between the TF and variational 
solutions as reported by the Eq. (1T7|) at A — > oo. Concerning the perturbation theory, the 
best accuracy, lees than 3%, is reached for A < 2. The results shown in the Figs. 3 and 
4 have a universal character and the comparison between the analytical methods provides 
universal criteria of their validity ranges. 



B. Attractive interaction 



Following the same trends as in the repulsive case, Fig. 5 shows the normalized order 
parameter for four negative values of A. In the figure, the dotted line represents the soliton 
solution ([5]). We observe that using the series (|24"|) the agreement is not so wide ranging 
as for the repulsive case . For values of A > —10, we obtain a better match between the 
Eq. ( 1241) and the numerical solution ( |33|) . Nevertheless, the agreement reached with the 
soliton solution ([5]) is remarkably good. In order to quantify the discrepancy between the 
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FIG. 7: (Color online) The same as Fig. [4] for the attractive interaction as a function of A. Dotted 
line represents the soliton chemical potential ([6]). Inset: Perturbation theory (|3ip (dot-dashed line) 
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FIG. 8: Effective Potential Veff for the GPE (see text). 

implemented analytical solutions and the numerical one, we evaluate the accumulated error 
( 136]) in terms of A. Figure 6 presents rji for all considered functions $j. Here, a dotted line is 
used for the soliton solution (j3J). We have estimated that the best result by using Eq. fpHj) is 
reached, for A > —10, while the exact soliton solution gives a better approach for A < —10, 
and in both cases we have an accumulated error rj < 0.005. The soliton variational solution 
®var yields a maximum error of T] ~ 0.36 at A « —4.6. The accumulated error using the 
perturbation wave function (dash-dot line) is also shown in the inset and as expected ?7 — > 
as A — > 0. 

Figure 7 depicts the calculated chemical potential /x for the variational calculation (Eqs. 
(JEJ) and (jHJ)), the solution following Eq. f l28l) . soliton solution (jBJ), and the numerical im- 
plementation for the GPE. The numerical procedure for the calculation of fi/hut was imple- 
mented in order to achieve a maximum uncertainty of 10~ 10 . According to the results of Fig. 
7, the system (126]) using 50x50 matrix reproduces quite well the chemical potential values in 
the interval A > —10 with an accuracy less than 1.2%, while for fi s , given by ([61) , the relative 
error tends to zero as A decreases. The best accuracy for the solution ([9l) is reached in the 
interval —3 < A < and fails for smaller values of A. In the inset, we show the calculated 
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chemical potential in the framework of a perturbation method, Eq. (1301) . and compared with 
the other four methods. Here, it can be seen the strong deviation of the soliton solution from 
the correct values for A > —2. However, no differences are observed between the numerical, 
perturbation method, and the calculations using f[2T)j) . Again as in the repulsive case, the 
results shown in Figs. 6 and 7 are of universal validity, giving an absolute estimation of 
the accuracy of each employed method as a function of a unique dimensionless parameter 
A. The present results teach us the way to get simple and exact analytical solutions for the 
GPE in the attractive interaction case in terms of A. Indeed, for A > — 10 using a small base 
(of the order of 50 oscillator wave functions) we obtain an accuracy of 10~ 8 for the chemical 
potential along with a minimum accumulated error, 77, of 0.005 for the order parameter. For 
smaller values of A the soliton solution and can be implemented as the exact solution 
of Eq. ([T]). At this point, it is necessary to analyze the convergence to the exact solution 
provided by the series (124"|) . In principle, as it was derived in Sec. II, the function (1241) is an 
exact representation of the order parameter <3> with a convergence at list in energy to the real 
order parameter $. The basis {ip n (x)} is a complete set for the Hilbert space defined by Eq. 
(CE]) independent of the sign of the non-linear interaction term. Nevertheless, the number of 
the harmonic oscillator wave functions needed to reach the necessary convergence to the real 
solution depends on the values and sign of A. The key point is to know when the series (124|) 
is really a good method for calculations and more efficient than the numerical ones. In our 
case, we selected 50 even functions (p n reaching an accuracy for the chemical potential less 
than 10 -8 in the range — 10 < A < 25. To get the same accuracy for the chemical potential 
in the attractive region with A < — 10, it is necessary to deal with matrixes (1281) of rank 
larger than 50x50. 

In order to clarify this peculiarity of the expansion ff24|) we define the effective potential 

Veff = ^mu 2 x 2 + A |$ nMm | 2 , 

where the order parameter $ has been substituted by the numerical solution $ num . Figure 8 
shows the potential Veff in units of ftw for both, the attractive and repulsive interactions. 
In the figure, we represented the exact calculation of fi for each considered value of A. It 
becomes clear that for the repulsive case, Veff resembles the harmonic oscillator potential 
(Fig. 8 a)) and the chemical potential falls within certain range of the harmonic oscillator 
eigenvalues. Hence, the complete set of harmonic wave function {{p n (x)} can reproduce well, 
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with an inexpensive computational effort, the mathematical properties of the GPE. In the 
case of attractive interaction, see Fig. 8 b), the situation changes drastically. Here, the 
effective potential becomes more localized as A decreases and for A — > — oo, Veff ~ S(x). 
The function Veff does not resemble the harmonic oscillator potential, thus the values of 
the chemical potential are far away from (n + |) eigenvalues. Although the basis {(p n (x)} 
is complete, the number of functions ip n (x) needed to describe the order parameter $ and 
chemical potential [i with certain accuracy should increase enormously as A decreases. This 
performance of the attractive interaction, determines that the Green function solution or 
equivalently the Galerkin or spectral method becomes computational expensive and the 
method is not adequate to describe the GPE for strong attractive interaction case, that is 
for A < -10. 

IV. CONCLUSIONS 

We have provided simple analytical forms to get explicit solutions for the GPE. The re- 
ported analytical techniques allow us to explore regions of positive and negative nonlinear 
interactions in condensates. We estimated the range of applicability of the perturbation 
theory, Thomas-Fermi approximation, soliton wave function, soliton variational calculation, 
and Green function solution (spectral method) through a universal interaction parameter 
A = X/IqTuv. The perturbation method is valid in the weak interaction limit, —2 lohw < A < 2 
lohu with an error for the chemical potential less than 1.5% while the TF approximation 
provides an error less than 3% if A ^ 10 l huj. The solution (|24p with solely 50 harmonic 
oscillator wave functions reproduces quite well the chemical potential /i with an accuracy of 
1.2% in the interval —10 l huj < A < 10 IqFujJ. We identified that the series fl24|) or the spec- 
tral method is not adequate and can be computational expensive for the attractive case if 
A < —10 lohu (see Figs. 6, 7 and 8). In this case, the bright soliton solutions ([5]) and (EI) rep- 
resent the better approach for the order parameter and the chemical potential respectively. 
The presented soliton limit is formally equivalent to the Thomas-Fermi one and becomes a 
powerful tool for condensates with strong attractive interaction. Also, we have introduced 
a soliton variational procedure valid for repulsive and attractive interactions which can be 
applied to the study of the dynamics of BEC or to model physical systems obeying the GPE. 
With the present results it is possible to have a short and comprehensive discussion on the 
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usefulness of different approaches for the mathematical and physical description of the BEC. 

We should note that the mathematical models here developed can be straightforward 
extended to the three-dimensional case, 9 two-dimensional "pancake-shaped",— or to the 
"cigar- shaped" BEC's^ 1 ^ 1 ^ and to study the dynamics of two component BEC systems.— 
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APPENDIX A: VARIATIONAL CALCULATION 

Inserting the wave function (jSj) in Eq. ([TJ) follows the Eq. (jTj), where a, /3, and 7 are 
numbers equal to: 




(Al) 



(A2) 



(A3) 



APPENDIX B: MATRIX ELEMENTS 



The fourth dimensional matrix T introduced in Eq. ff27j) is defined as 



1 



plmn 




(Bl) 
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The matrix elements T p i mn have the followings properties: 

i) Tpi mn = if n + m + I + p =odd number. 

ii) T p i mn is invariant under the permutation of the quantum numbers n, m, I, and p, i.e. 

Tplnm T[p mn Tp m i n ... 

iii) For m = we find^i 

_ 2'- 1 r( a -or( a - P )r( a -n) 

^+H\p\n\ ' 1 ' 

where T(z) is the gamma function and 2s = n + l+ p+ l. 

iv) The following relations hold between two successive matrix elements T p i n0 : 

T - (s-l-l)(3-n-l) 

J-plnO — ; — 1 : J p-2/n0) V a6 ) 

(s-p)VP(P- 1 ) 



or 



TpJnO — ~ 7F= ~^p-l£-lnO> (B4) 



with 



v 



2tt 

) For the most general case we have the expressio 



^oooo — 7= • (B5) 
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- 1 1 M ~ rn ~P'2 M ~ \ 



- p,l,n,m 



X 



ny/2 n + m + l +Pn\m\l\p\ 
T{M -l + \)T{M -n+\) 
T{M-n-l + \) 



x 



F . ~m, -p, -M + n + l + l;^ 
-M + Z + |, -M + n+i; 



(B6) 



. «i, «2, Ct3; i . ■>] ■ 

where 3F2 1 is the generalized hypergeometric series^ and 2M = p + I + 



Pi, A 



2. 



m + n. 

vi) The matrix element T p i nm satisfies the recurrence relation 
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-Lplnm \ J-pln+lm— 1 

v. m 

— Tpl-lnm-1 — \ —Tp-Unm-l- (B7) 

m \ m 

These mathematical properties allow to evaluate the tensor T in a straightforward way and 
in consequence to solve Eq. ( 1281) for the eigenvalues /i and eigenvector C very efficiently 
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